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Due to the time scale problem, rare events are not accessible by straight forward molecular 
dynamics. The presence of multiple reaction channels complicates the problem even further. The 
feasibility of the standard free energy based methods relies strongly on the success in finding a proper 
reaction coordinate. This can be very difficult task in high-dimensional complex systems and even 
more if several distinct reaction channels exist. Moreover, even if a proper reaction coordinate can 
be found, ergodic sampling will be a challenge. In this article, we discuss the recent advancements 
of path sampling methods to tackle this problem. We argue why the path sampling methods, via 
the transition interface sampling technique, is less sensitive to the choice of reaction coordinate. 
Moreover, we review a new algorithm, parallel path swapping, that can dramatically improve the 
ergodic sampling of trajectories for the multiple reaction channel systems. 
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I. INTRODUCTION 

Path sampling has shown to be efficient tool to study 
rare events that are not accessible by straight forward 
molecular dynamics (MD). The principal idea behind the 
method [T] is to reduce the superfluent exploration of the 
stable states and to generate trajectories that have a high 
chance to become reactive. Via the transition interface 
sampling (TIS) method 0, H, 0, sets of trajectories 
can be generated that progressively climb across the bar- 
rier towards the product state. Besides the rate of the 
reaction, the TIS path ensembles allow to analyze the 
mechanism. Each TIS simulation in the series give a cor- 
rect distribution of pathways that reach a certain level 
towards reactivity. Some of them will fail to reach a next 
level and some of them will successfully make a step fur- 
ther. Therefore, these path ensembles yield a treasure of 
information for analyzing complex reaction mechanisms 
in a very intuitive way. The main TIS equation relates 
the reaction rate with a flux through a surface closeby 
the reactant state times the overall crossing probability 



kAB = /a7'a(As|Aa)- 



(1) 



The flux fA is simply the number of crossings with sur- 
face \a per unit time. 'Pa{^b\^a) is the probability that 
whenever the surface Xa is crossed, the system will go 
all the way over the barrier to cross the surface \b at 
the other side. The surfaces are generally defined by a 
parameterization using a single parameter called reaction 
coordinate (RC) or order parameter. The RC can be any 
non-linear function of all particle positions and velocities 
in the system js^. The surface (or interface) A is then 
simply defined as the collection of phase points for which 
this function is exactly A. Typical examples of RCs are 
distances of bonds that have to broken or formed, the 
size of the largest solid cluster in nucleation studies, or 
the number of native bonds for protein folding. By con- 
vention, we assume that \a < A^. It is important to 
realize that the reaction rate k does not depend on the 
choice of the function RC or on the values A^i and 



as long as they obey some very basic principles. That is, 
each trajectory from A^ to As should be a true reactive 
event of the reaction of interest. In addition, the chance, 
that the system quickly returns to Xa after As is crossed, 
should be of the same order as that of an independent 
reverse reaction. 

/a can be computed by straight-forward MD. 
'Pa (As I A^) is very small and therefore difficult to com- 
pute. However, if we define a set of interfaces in between 
with Aq = Ayi,A„ = Xb, and < Ai_i, we can invoke 
following exact factorization: 

Va{Xb\Xa) = VA{Xn\Xo) ^Y[VA{X, + l\Xi). (2) 

1=0 

Here, ^^^(Ai+i | AJ is a history dependent conditional 
crossing probability that is much higher than Pa(Ab |Aa) 
and, therefore, much easier to compute. It equals the 
probability that, given the system is about to cross Ai 
for the first time since its last crossing with Xa, it will 
cross Ai-i-i as well before A^. Computing this factor can 
be done by generating a representative set of trajecto- 
ries that start at Xa, cross Ai, and end at either A^ or 
Ai+i. The fraction that ends at A^+i equals VAiXi+i\Xi). 
Recent TIS simulations employ a somewhat different en- 
semble. They consider all possible pathways that start 
at Xa, end at either Xa or As, and have at least one 
crossing with A^. 7'yi(Ai+i |Ai) is then the fraction that 
crosses A^+i. Generating this set is only slightly more 
expensive as most trajectories will end at Xa anyway. 
However, this approach makes it much easier to initialize 
the TIS parameters, such as the interface positions, while 
running one path simulation after the other in the sim- 
ulation series. Moreover, the implementation of the new 
path swapping technique is much more straightforward 
when these ensembles are used. Of course, the crucial 
point is how to generate these trajectories. In principal 
we could simply run a MD trajectory and cut out the 
pieces of interest, but then our efficiency would be as 
bad as MD. Luckily, Dellago et al. developed efficient 
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Monte Carlo moves to generate these trajectories within 
the context of transition path samphng (TPS) The 
algorithm only requires to provide a single trajectory that 
satisfies the conditions. The principal MC move, called 
shooting, picks a random timeslice of the old trajectory, 
changes the velocities by a small amount, and integrates 
the equations of motion backward and forward in time to 
obtain a new trajectory. This trajectory can be accepted 
if it also satisfies the opposed conditions. The MC move 
in TIS is a slight adaptation of this shooting move to al- 
low flexible path lengths (for a graphical illustration of 
the TIS shooting algorithm see [7|). This flexible path 
length is a strong improvement upon the original TPS 
algorithm. The trajectories in TIS will on average be 
shorter than in TPS. In addition, TIS docs not introduce 
any systematic error as the small fraction of long trajec- 
tories is also taken into account. These will usually lie 
outside the range of path lengths one would consider in 
TPS. Other improvements are that the TIS formulation 
of the reaction rate is less sensitive to recrossings and 
that TPS shifting moves have become redundant in the 
TIS algorithm. 

The TIS algorithm has been successfully applied to 
several systems ranging from protein folding [8], nucle- 
ation fa], rnicelle fusion and fission [13], and DNA dcnat- 
uration [ll|. Moreover, the TIS theory has initiated the 
development of several new path sampling methods such 
as the partial path TIS (PPTIS) [l2| and forward flux 
sampling (FFS) 13]. PPTIS was especially designed for 
diffusive barrier crossings. It uses a history-loss assump- 
tion to reduce the path length even further. FFS is exact, 
like TIS, but uses another type of MC move in which 
one integrates the equations of motion only forward in 
time . The advantages of FFS is that it works for both 
equilibrium and out of equilibrium systems and that the 
dynamics do not necessarily have to be reversible. The 
disadvantage is that FFS introduces much stronger corre- 
lations between data points. This has a negative effect on 
the efficiency of the method, which can become extreme 
if the dynamics bear a strong deterministic character or 
when the RC is not weh chosen [11]. Both PPTIS and 
FFS are based on Eqs. Q and ©. 

In this article we will discuss the computational chal- 
lenging problem of multiple reaction channels. This prob- 
lem is quite generic for any transition that involves many 
important degrees of freedom and it poses several diffi- 
culties. First of all it is extremely difficult to derive good 
RCs (like committor surfaces [l^, [1^1 ) for such a system 
that can be efficiently be used in a free energy based 
approach. Therefore, one would either need a system- 
atic approach to obtain feasible RCs or a method that is 
much less sensitive to the choice of RC than the standard 
methods. Secondly, in order to explore one channel after 
the other, one needs a MC sampling that is highly nonlo- 
cal. Finally, the method should be very flexible and able 
to let the nonlocal MC moves coincide with a high rate of 
acceptance. We will show that TIS encompasses both the 
required insensitivity to the RC and the nonlocal char- 



acter of shooting moves. Moreover, using the additional 
new technique of parallel path swapping (PPS), we can 
also fulflU the last requirement upto a very satisfacto- 
rily level. This results in a very efficient approach to 
treat multiple reaction channel systems. This article is 
organized as follows. In Sec. |TT1 we discuss the issue of 
RC sensitivity for TIS and related methods. In Sec IIIII 
we introduce the new PPS technique and review recent 
results, based on this method, for the DNA denatura- 
tion [ll|. We end with conclusions in Sec. IIVI 



II. REACTION COORDINATE DEPENDENCE 

The problem of finding suitable RCs in high- 
dimensional complex systems has been an outstanding 
problem for several decades. The free energy based ap- 
proaches, pioneered by Wigner [13], Eyring[3l, and Keck 
[l^ and further developed by Bennett jl^] , Chandler [2l[ , 
rely on a two-step approach. First the free energy as func- 
tion of the proposed RC needs to be calculated using um- 
brella sampling (US) [22i] or thermodynamic integration 
(TI) [2^ techniques. From this, the free energy barrier 
can be derived and the transition state theory (TST) ap- 
proximation for the rate. For complex systems, the TST 
approximation needs to be corrected for fast recrossing 
events by a dynamical factor k (with < k < 1), called 
the transmission coefficient. If k is not too low, it can 
be accurately determined by releasing dynamical trajec- 
tories forward and backward in time from the top of the 
free energy barrier, n is then related to the fraction of 
trajectories that connect reactant state A with product 
state B [3l| . In this way, one can determine k exactly, 
independently to the choice of RC, and one obtains a 
set of reactive trajectories that can be subject of further 
investigations. However, if the RC fails to capture the 
actual mechanism, the vast majority of trajectories are 
either of the type A ^ A or B ^ B and k becomes im- 
measurably small. Illustrative for this effect are chemical 
reactions in solution which require considerable solvent 
rearrangements. The top of the free energy barrier, de- 
fined by a RC that does not include the solvent degrees of 
freedom, will either correspond to the situation where the 
solvent can easily incorporate the reactant or the product 
species. As result, both forward and backward dynam- 
ics will rapidly collapse to the same state. In addition, 
the evaluation of the free energy barrier itself becomes 
problematic. Considerable hysteresis will occur in US 
or TI when the system is dragged over the barrier from 
reactant to product state and back. 

This RC problem was the main driving force to de- 
velop alternative methods like TPS [ij. Although it is 
quite evident that TPS can generate reactive trajectories 
quite efficiently using a simple type of RC, it is not so 
clear if the same is true for the actual reaction rate eval- 
uation. Only recently, Ref. [ij] really proved, by analyt- 
ical results of a simple 2D model, that the TIS efficiency 
of the reaction rate calculation is much less sensitive to 
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an improper choice of the RC than free energy based 
methods. There are two reasons for this. First of all, 
as paths are global identities (the weight of a path does 
not depend on a single point), the sampling of paths in- 
stead of phase points can eliminate hysteresis effects [ill ■ 
Secondly, whereas free energy based methods apply im- 
portant sampling approaches to configuration space only, 
hoping that the dynamical factor will be not too small, 
path sampling applies importance sampling techniques 
on the dynamical factor directly. The RC insensitivity is 
not generic for any type of path sampling method. It is 
valid for the original TPS method as well, but not for e.g. 
PPTIS or FFS JJj. We can also give a third reason which 
was not important for the 2D example studied in [3] ■ As 
the TIS/TPS shooting move is nonlocal, it should be able 
to circumvent barriers that are orthogonal to the RC. In- 
deed, a TPS water trimer study [2J] revealed that the 
shooting algorithm was capable of finding two reaction 
mechanisms across different saddle points separated by a 
barrier higher than the total energy of the NVE simula- 
tion. An approach to enhance this effect is explained in 
the next section. 

The different dependence on the RC of TIS and FFS, 
although based on the same equations, can be explained 
by the argument that Eq. ([2]) has two possible inter- 
pretations. The conditional crossing probabilities can 
be viewed as a kind of physical (non-Markovian) hop- 
ping probabilities to go from one interface to the other 
like climbing up a ladder. This interpretation is very 
close to the FFS implementation of Eq. ^ . The al- 
ternative interpretation is that Eq. ^ simply corre- 
sponds to a multiplication of ratios between the num- 
ber of paths contained by different sets: VAi^Bl^A) = 

#paths g[n+] _ #paths g[l+] #paths g[2 + ] #paths g[3+] 
#patlisg[0+] ~ #pathsg[0+] #pathsg[l + ] ^ #pathsg[2+] ^ 



Product State B 



[0 ] Brute Force 
1000,000 trajectories 



TIS Simulations 100 trajectories each 



:^paths g[ 



where 



defines the collection of 



■ ■ ■ ^ #pathsg[(n-l) + l 

paths that start at and have at least one crossing with 
Xi before revisiting Aa or ending at As. From this, it is 
directly clear that if we replace e.g. [3+] in this factor- 
ization by an arbitrary different set of trajectories, this 
still does not change the validity of the equation. This 
has an important implication. The final trajectory set 
[n~^] can be fundamentally different than [3+]. They do 
not necessarily have to resemble up to A3. 

In Fig. [1] we give two simple examples for seven in- 
terfaces {Ao, . . . , Ae}. Shown here are path survival di- 
agrams. In Fig. [l]-(top), we assume that only 1 out of 
the 10^ trajectories that initially start at A^ will reach 
A_B and at each interface only 10% survives reaching the 
next one. Hence, if we would simply run straightforward 
MD trajectories, we would need at least 10^ trajecto- 
ries to hope for a single reactive event. Now, if we only 
run 100 trajectories per TIS (or FFS) ensemble, we ex- 
pect that each time 10 trajectories will reach the next 
level. Our final result ^^(AaIAb) = ]Xi=i 'PAiKlXi-i) = 
(10/100)^ = 1-10~^ is exactly identical to a perfect brute 
force calculation but using only 600 trajectories instead 
of 10^ and having 10 fully reactive trajectories instead 
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FIG. 1: (color online) Path survival graphs. Top: a possi- 
ble case for a one channel system. Interfaces are positioned 
such that one out of ten reaches the next interface. The left- 
hand side is the situation when trajectories are generated by 
brute force using one million trajectories. The right-hand side 
shows the results (assuming a perfect sampling) for the TIS 
simulation series, each simulation consisting of only one hun- 
dred trajectories. Bottom: a possible case for a two channel 
system. Again a brute force simulation of one million trajec- 
tories is compared with a TIS simulation series of one hundred 
trajectories each. One channel (red) has initially a higher sur- 
vival rate than the green one, but turns into a dead end. The 
paths gathered by TIS for the [2"*"] ensemble do not contain a 
single green path. Still, the product of the obtained crossing 
probabilities gives the right result. 



of only 1 in the end. This directly shows the orders 
of magnitude improvement of the path sampling simu- 
lations compared to MD. In a bit more elaborate cal- 
culation, one can show that the effective computational 
cost scales exponentially with the barrier height in brute 
force MD and only quadratically in TIS^ This is similar 
to US using rectangular bias windows 14). In Fig. [l]- 
(bottom) , we show a bit more complicated example with 
two channels. One is initially favorable but turns into a 
dead end, while the other one survives. The right-hand 
side shows the most likely outcome for the TIS simula- 
tion series. It is important to note that at some point, 
at the [2+] simulation, the set of trajectories generated 
by TIS does not contain a single green path. Still the 
final result ^^(AelAo) = (10/100) ■ (20/100) • (10/100) • 
(6/100) • (8/100) ■ (10/100) = 0.96 • IQ-^ is correct. This 
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FIG. 2: (color online) Illustration of the multiple channel bar- 
rier. Contour plot is shown below. Trajectories a and c are 
two possible trajectories in the [1^] ensemble situated in two 
different channels. In order to find trajectory c from trajec- 
tory a, we would either need to shoot (with some luck) from 
inside the reactant well, or to generate the bridging trajectory 
b that is high in energy and thus likely rejected. Trajectories 
d and e are trajectories in the [3^] ensemble. They are much 
higher in energy and can easily move from one channel to the 
other. 

is a very important point to realize as it shows once more 
the flexibihty of the TIS method. FFS which propagates 
trajectories only forward in time will not be able to re- 
capture the green trajectories once they are lost. Hence 
FFS requires to sample much more trajectories at the 
initial interface ensembles or to device a RC where this 
problem will not occur. Of course, in order to recapture 
the successful green pathway it is of eminent importance 
that the TIS shooting move will be able to make this 
transition. In the next section we show how PPS can be 
a very effective tool to enhance these transitions. 

III. PARALLEL PATH SWAPPING 

In Fig. [5] we depict the problems that can occur when 
the standard shooting move is applied for a multiple re- 
action channel system. The trajectories a and c of the 
[1+] ensemble are situated in different channels of the 
potential energy barrier. It is not very likely that the 
shooting move will be able to connect these trajectories. 
It would require to generate a bridging trajectory b, but 
such a trajectory is very high in energy and will likely 
be rejected. Alternatively, in order to make use of the 
non-locality of the shooting move and its ability to avoid 
transversal barriers, we might want to shoot from inside 



the reactant well. However, such a point is not part of the 
paths of this ensemble that start at A^. We might shift 
the Xa interface more into the reactant well with the ex- 
pense that all trajectories become longer and thus more 
expensive. However, even then we must be very lucky 
because if we fail to cross Ai, this trajectory must be re- 
jected and we remain in the same channel. On the other 
hand, if we consider the trajectories of the [S"*"] ensemble, 
these are much higher in energy and can easily move from 
one channel to the other. It would be very useful if we 
could somehow make advantage of the high energy paths 
of the [S+J ensemble and the non-local shooting moves 
from inside the reactant well. 

A very successful method in standard MC is the par- 
allel tempering or replica exchange method [25| . In this 
method one performs several simulations in parallel at 
different temperatures. Then, with a certain frequency 
and acceptance probability the configurations at one tem- 
perature simulation is being swapped with a lower tem- 
perature simulation. The high temperatures will easily 
explore the rough free energy surface as they have a much 
lower probability to get trapped. The low temperature 
simulations will benefit from the exchange of information 
as they will be able to hop from one potential basin to 
another without having to cross the intermediate barriers 
physically. The combination of path sampling and par- 
allel tempering has been used before [26], but the main 
disadvantage of this approach is that one needs to per- 
form expensive additional simulations whose information 
is useless if one is interested in the reaction rate at one 
temperature alone. In addition, parallel tempering will 
not help to circumvent entropic barriers using the non- 
locality of the shooting move. However, if one realizes 
that the TIS method already consists of a series of sim- 
ulations, it makes sense to introduce a swapping move 
between these [llj . 

In order to create the maximum flexibility of swapping 
moves at all levels, it is convenient to replace the initial 
MD simulation for the flux calculation /a by another 
type of path ensemble [0~]. These consists of all the 
paths that start at Xa, then go initially towards the neg- 
ative direction, and finally end again at Xa- The initial 
flux can than be obtained from the average path lengths 
in the [0~] and [O"*"] ensemble: 

/.= (0 + 0)", (3) 

where (ip^tlj), (^path) ^^'^ average path lengths in 
these ensembles. In this series of path ensembles, 
{[0~], [0+], . . . ,[{n — 1)"'"]}, the swapping move be- 
comes extremely effective. Note that the swapping moves 
do not require any force calculations except for the swap- 
ping between [0~] and [0+] (see Fig. [3]). Here, the last 
timestep of the old path in the [0~] ensemble is used 
as initial point to generate a new trajectory in [0"*"] by 
integrating the equation of motion forward in time. Con- 
versely, the initial point of the old path in [0+] is followed 
backward in time to generate a path in [0~]. The TIS 
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FIG. 3: Illustration of the swapping moves [0~] <~* [0^] and 
[1"*"] <-» [2+]. The alternative swapping possibility [0^] ^ [l"*"] 
would have yielded a rejection as the [0"*"] path does not cross 
Ai and is therefore not a valid path for the [1^] ensemble. 



algorithm is then as follows. At each step it is decided by 
an equal probability whether a series of shooting or swap- 
ping moves will be performed. In the first case, all simu- 
lations will be updated sequentially by a shooting move. 
In the second case, again an equal probability will decide 
whether the swaps [0"] ^ [0+], [1+] ^ [2+],... or the 
swaps [1+] ^ [2+], [3+] ^ [4+], . . . are performed [ll|- 

The effectiveness of this algorithm was illustrated 
in [llj for the denaturation transition of DNA using the 
mesoscopic Peyrard-Bishop-Dauxois (PBD) model [13]. 
In this model, the DNA molecule is represented by a 
sequence of one dimensional particles representing the 
relative base-pair separations from the groundstate posi- 
tions. Each particle is positioned in an external Morse 
potential describing the interaction of base-pairs of op- 
posite strands. In addition, a first-neighbor anharmonic 
spring potential is used for the stacking interaction be- 
tween bases of the same strand. The width and depth 
of the Morse potential are adapted to describe the weak 
AT or the strong GC interaction. Due to thermal fluctua- 
tions the hydrogen bonds between base-pairs of opposite 
strand can break, which corresponds to a particle moving 
on the plateau of the Morse potential. However, if the 
neighboring particles are still in the closed state, this par- 
ticle will be rapidly pulled back into the stack. The fully 
denaturated state is achieved when all base-pairs move on 
Morse plateau after which the two DNA strands have no 
interaction anymore and can move to infinite distances. 
This event is very rare for the larger molecules and has 
a very complex dynamics as it can proceed via different 
path ways. The DNA molecule might initially open up at 
one end and propagate the opening through the molecule. 
Alternatively, a bubble in the middle might appear that 
continues to grow in both directions. Henceforth, the 
denaturation process is a typical example of a multiple 
reaction channel system and an accurate evaluation of 
the rate is quite a challenge for the larger molecules. 



In Ref. the dynamics of a 20 base-pair DNA 
molecule of AT bases was investigated by TIS with and 
without swapping. As RC, the base-pair separation of 
the base-pair with the smallest distance was used. In to- 
tal eight interfaces {Ao, . . . , Ay} were defined (for more 
details see Ref. [HI). The calculated rates rates were 
in good agreement 0.0492 ± 0.0062 and 0.0524 ± 0.0025 
ns~^ for standard TIS and TIS with path swapping re- 
spectively. The latter has a significant lower error despite 
a much shorter simulation. A very accurate integration 
method (28j for quasi one-dimensional systems confirmed 
the path swapping result within a 0.6 % uncertainty. 

A very instructive approach to quantify the efficiency 
of the individual simulations and the simulation in total 
is given in [l^l by the introduction of so-called efficiency 
times Tcff . These are defined as the number of force calcu- 
lations required to obtain a statistical error equal to 1 in 
a given simulation. For the TIS simulations [O+J, . . . 
these can be expressed as fl^ : 



1 -Pz 

Pi 



(4) 



Here, = Va{Xi+i\Xi) and L, = (tf^'O/At with At the 
MD timestep. These are in principle independent to the 
simulation method, is the ratio between the average 
cost of the simulation cycle and Li. Mi is the effective 
correlation between trajectories. 

Fig. m shows the five parameters in Eq. ^ for the 
seven TIS path ensembles using standard TIS and path 
swapping. The results for pi and Li are the same as ex- 
pected. The values for ^i show that the average cost per 
cycle is reduced by a factor 1.5 for the [0"''] and by a fac- 
tor 2 for the [i"*"],? > 0, simulations. This is due to the 
swapping moves that do not require any force calcula- 
tions except for [0~] <-> [0^]. However, more importantly 
is the dramatic reduction of A/i . The large values of Mi 
in the standard TIS simulations are directly related to 
the problems of ergodic sampling when the system gets 
stuck for a long time in one specific reaction channel. 
Most spectacular is the reduction in the [0+] ensemble 
by more than two orders of magnitude. The reduction in 

both and Mi is reflected in t^^ ' . The computational 
cost is reduced at all levels by at least a factor of two, 
but much more for ensembles [0+], . . . , [4+] that consists 
of shorter paths with lower energy. When all results are 
taken together, it can be concluded that the path swap- 
ping technique resulted in a gain of efficiency by more 
than a factor 20. 



IV. CONCLUSIONS 

We have discussed the ability of path sampling to study 
transitions that proceed via multiple reaction channels. 
This is a common phenomenon for any complex reaction 
mechanism that involves many degrees of freedom. The 
analysis of these processes are a huge challenge from a 
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FIG. 4: (color online) Efficiency analysis for the TIS inter- 
face simulations with (green) and without (red) path swap- 
ping. The crossing probability pi, the average discrete path 
length Li, the ratio and the effective correlation A/i (note 
the change of scale in the y-axis) is shown in the first 4 pan- 
els. The last panel shows the effective computational cost per 
interface simulation. For the actual data see Ref. JT| 



computational point of view. Besides all the difficulties of 
rare event simulations, many additional issues come into 
play. The development of a feasible RC, that can be used 
in a standard free energy based method, is notoriously 
difficult for high dimensional complex systems. Finding 
a RC that can well describe several reaction mechanisms 
simultaneously is even more difficult than that. 

Alternatively, one could rely on a method for which the 
correct RC does not play such a crucial role. We have 
given several justifications why this is the case for the TIS 
and TPS path sampling methods. For a 2D model system 
of a slanted barrier, the advantageous scaling of TIS com- 
pared to standard methods was proven rigorously when 
the effective computational cost as function of the 'qual- 



ity of the RC was calculated analytically [l^. In sec. HIl 
we gave some additional arguments based on a fictitious 
two channel problem displayed by the path-survival dia- 
gram of Fig. [1] The difficulty of this example is that one 
channel is initially favorable but finally turns into a dead 
end. If we assume a perfect ergodic sampling, the TIS 
simulations are still able to return the right result using 
a minimal amount of trajectories. FFS, in which trajec- 
tories are propagated only forward in time, would not be 
able to achieve the same. It can not change the history of 
the paths and will therefore get stuck in the unfavorable 
channel. Apart from this, to be able to sample between 
several distinct reaction mechanisms, one needs to com- 
bine nonlocal MC moves with a high rate of acceptance. 
The shooting move has shown to have the required nonlo- 
cal characteristics [l^l , but in order to sample efficiently, 
transitions between the different channels need to occur 
at a much higher frequency than in the standard shooting 
algorithm. 

Sec. lIIII shows a very promising approach to accomplish 
this based on the ideas of parallel tempering [2^ . Instead 
of performing several path sampling simulations at differ- 
ent temperatures |26|, the swapping occurs between the 
different interface ensembles. The initial MD simulation 
to compute the flux is replaced by an 'internal' path en- 
semble which allows to swap trajectories at all levels in 
the system. Compared to standard parallel tempering, 
this gives the advantage that there is no need to extend 
the number of TIS simulations. It also provides a way 
to overcome entropic barriers for which parallel temper- 
ing would not help. We reviewed the results of Ref. [ll| 
where this method was applied on the DNA denaturation 
transition using the mesoscopic PBD [2^ model. These 
results showed that the PPS technique improved the TIS 
efficiency by more than a factor 20. 

Still, there are many issues that need to be studied and 
more improvements might be possible. An exact estimate 
of the total error becomes more complicated as the stan- 
dard error propagation rules assume that the different 
simulations are independent. We think that the presence 
of covariant terms will not have a huge effect, but this 
will be investigated explicitly in the near future. In p3 | 
we derived relations for how to divide a total fixed simu- 
lation time over a simulation series to obtain the lowest 
possible overall statistical error. The optimum was found 
when each simulation was given a simulation time pro- 
portional to oc '\J ' . This result is generic for any 
type of method for which the final quantity is obtained 
by a product outcomes from a series of independent sim- 
ulations. Finding the optimal way to divide the total 
simulation time in a PPS algorithm is much more com- 
plicated. The simulations are no longer independent and 
each simulation has a double function. Besides the eval- 
uation of the intrinsic property, their function is also to 
assist the ergodic sampling of the other simulations. So 
far, the PPS simulations were performed using the same 
number of trajectories for each ensemble. We are now de- 
veloping an improved PPS-TIS method that iteratively 
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adapts the number of trajectories per ensemble to the op- 
timal ratios [29|| . We think that TIS in combination with 
PPS can become an important standard method for the 
sampling of rare events in complex systems and multiple 
reaction channels in particular. 
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